load c:\manoj\projects\plasma\HUA_PIU HUA_PIU_fday HUAPIU magnetic_frac;
data11 = load('c:\manoj\projects\plasma\STRICT_CRT_PLP_EOP_JT_09-17LT_268-298LON_KP.txt');
load c:\manoj\projects\plasma\Julia_PLP.mat;
bin_radius = 30/(60*24);%<------- NORMAL
%bin_radius = 60/(60*24);
kp_limit = 3;
lt_step =  1;
%for delay_index = 1:13,
delay_index = 1;
delay_sel = [13,14,15,16,17,18,19,20,21,22,23,24];
Kp_delay = delay_sel(delay_index);%Valid columns [10,13-24];
kp_delay_text = ['00';'01';'02';'03';'04';'05';'06';'07';'08';'09';'22';'MN'];
lag_step = 30;%< NORMAL 30
min_data = 5; %Minimum number of data pairs to be considered for correlation
LAG = [30:lag_step:390]/(60*24);
%LAG = [30:lag_step:240]/(60*24);%<------- NORMAL
LT = [10:lt_step:17]; %<------- NORMAL
%LT = [9:lt_step:17];
robft = 1;
plot_option = 1;%[1 scatter, 2 image 3 curves];
ED_matrix = ones([length(LT)-1,length(LAG)])*NaN;
LV_matrix = ones([length(LT)-1,length(LAG)])*NaN;
RT_matrix = ones([length(LT)-1,length(LAG)])*NaN;
PT_matrix = ones([length(LT)-1,length(LAG)])*NaN;
ED_error = ones([length(LT)-1,length(LAG)])*NaN;
LV_error = ones([length(LT)-1,length(LAG)])*NaN;
RT_error = ones([length(LT)-1,length(LAG)])*NaN;
PT_error = ones([length(LT)-1,length(LAG)])*NaN;

for jj = 1:length(LAG),
for kk = 1:length(LT)-1,
    
%      L = ~isnan(PLP(:,27)) & ~isnan(PLP(:,18))  &  D_KP(:,Kp_delay) <=kp_limit & PLP(:,7) >= LT(kk) & PLP(:,7) <= LT(kk+1)...
%          & PLP(:,12) == 1;% & PLP(:,20) >= 0 & PLP(:,29) >= 0;% & PLP(:,8) > -85 & PLP(:,8) < -70;
     L = ~isnan(PLP(:,27)) & ~isnan(PLP(:,18))  &  D_KP(:,Kp_delay) <=kp_limit & PLP(:,7) >= LT(kk) & PLP(:,7) <= LT(kk+1)...
         & data11(:,12) == 1;% & PLP(:,20) >= 0 & PLP(:,29) >= 0;% & PLP(:,8) > -85 & PLP(:,8) < -70;

         IndData = PLP(L,:);
         counter = 0;
         
         for i = 1:size(IndData,1),
         fday_champ = tmjd(IndData(i,1),IndData(i,2),IndData(i,3),IndData(i,4),IndData(i,5),IndData(i,6),1);   
         iday = find(HUA_PIU_fday==floor(fday_champ));
         if ~isempty(iday),
            fday_champ_frac =  fday_champ - floor(fday_champ) - LAG(1);
            LLL = isnan(HUAPIU(iday,:));kkk = find(LLL == 0);
            if length(kkk) > 3,
             if fday_champ_frac > magnetic_frac(kkk(1)) & fday_champ_frac < magnetic_frac(kkk(end)),
                 counter=counter+1;
                [trash left_bin] =  min(abs(magnetic_frac(kkk(1):kkk(end)) - (fday_champ_frac-bin_radius)));
                [trash right_bin] = min(abs(magnetic_frac(kkk(1):kkk(end)) - (fday_champ_frac+bin_radius)));
                left_bin  = kkk(1) + left_bin - 1;
                right_bin = kkk(1) + right_bin -1;
                HUA_data(counter) = nanmean(HUAPIU(iday,left_bin:right_bin));%<------- NORMAL
                %HUA_data(counter) = nanmean(HUAPIU(iday,860:1060)); % Daily
              
                Plp_ed(counter) = (IndData(i,18) + IndData(i,27))/2;
                Plp_lv(counter) = (IndData(i,17) + IndData(i,26))/2;
                Plp_ratio(counter) = (IndData(i,20) + IndData(i,29))/2;
                Plp_pt(counter) = Plp_ed(counter) -IndData(i,36);
                Plp_altitude(counter) = IndData(i,9);
             end;
            end;
         end;
         end;
        if counter >= min_data,
         LP = isnan(HUA_data);
        HUA_data(LP) = [];Plp_ed(LP) = [];Plp_lv(LP) = [];Plp_ratio(LP) = [];
        Plp_pt(LP) = [];
        counter = counter - sum(LP);
        end;
        if robft == 1 & counter >= min_data,
                [HUA_data,Plp_pt,Plp_ratio,Plp_lv,Plp_ed] = ...
                    rob_regress(HUA_data,Plp_pt,Plp_ratio,Plp_lv,Plp_ed);
           counter = length(HUA_data);
         end;
         lt_count(kk,jj) = counter;
         if counter >= min_data,
             [cof_ed,err_ed] = corrcoef(HUA_data,Plp_ed); 
             [cof_lv,err_lv] = corrcoef(HUA_data,Plp_lv); 
             [cof_rt,err_rt] = corrcoef(HUA_data,Plp_ratio); 
             [cof_pt,err_pt] = corrcoef(HUA_data,Plp_pt); 
             ED_matrix(kk,jj) = cof_ed(2,1);
             LV_matrix(kk,jj) = cof_lv(2,1);
             ED_error(kk,jj) = err_ed(2,1);
             RT_matrix(kk,jj) = cof_rt(2,1);
             RT_error(kk,jj) = err_rt(2,1);
             PT_matrix(kk,jj) = cof_pt(2,1);
             PT_error(kk,jj) = err_pt(2,1);
             LV_error(kk,jj) = err_lv(2,1);
             
             if plot_option == 1,
             subplot(length(LAG),length(LT)-1,(jj-1)*(length(LT)-1)+kk);
             %set(gca,'FontSize',16);
             hold on;
             plot(HUA_data,Plp_pt/1e6,'r.','MarkerSize',15);
             if jj ~=length(LAG), set(gca,'XTick',[]);end;
             if kk ~=1 ,set(gca,'YTick',[]);end;
             %grid on;
             axis([-30 150 0 2]);
             text(-20,2.5,sprintf('%d',LT(kk)));
             text(30,2.5,sprintf('%5d',LAG(jj)*(60*24)));
             text(150,2.5,sprintf('%4.2f',cof_pt(2,1)),'color','r');
             %xlabel(sprintf('%5d-%5d',LT(kk),LAG(jj)*(60*24)));
             end;
        end;
%         fprintf('\n');
         clear HUA_data Plp_ed Plp_lv Plp_ratio Plp_pt;
      
end;
end;

for i = 1:length(LT)-1,
SS(i,:) = sprintf('%2d-%2d',LT(i),LT(i+1));
end;

if plot_option == 3,
color_index = ['k*-'; 'b*-' ;'r*-'; 'g*-'; 'c*-' ;'y.-' ;'m.-'; 'r>-' ;'g.-'; 'c.-'];
AA = LAG*(24*60);
% color_map = colormap('jet');
% color_index = color_map([1:3:64],:);
figure(1);
subplot(331);
for i = 1:length(LT)-1,
hold on;
errorbar(AA,PT_matrix(i,:)',PT_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

title('C-T - Julia correlation');
text(-42,1.0,['Kp-' sprintf('%d',kp_limit)]);
text(-42,0.8,['Delay-' kp_delay_text(delay_index,:)]);
subplot(334);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(337);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(PT_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);


subplot(332);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),RT_matrix(i,:)',RT_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('CTR - Julia correlation');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

subplot(335);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

subplot(338);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(RT_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);


subplot(333);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),LV_matrix(i,:)',LV_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('L-Value - Julia correlation');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

subplot(336);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

subplot(339);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(LV_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');
set(gca,'XTick',AA);
set(gca,'XTickLabel',AA);

end;

if plot_option == 2, %image
    figure(2);
    AA = LAG*24*60;
    BB = [LT(1:end-1)+diff(LT)/2];
    subplot(2,2,1);
    imagesc(BB,AA,RT_matrix');
    set(gca,'YTick',AA);
    set(gca,'YTickLabel',AA);
for i = 1:length(AA),
for j = 1:length(BB),
text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
end;
end;
    
    
    caxis([0,1]);
    title('CTR');
    colorbar;
  %  set(gca,'FontSize',16)
    subplot(2,2,2);
    imagesc(BB,AA,PT_matrix');
        set(gca,'YTick',AA);
    set(gca,'YTickLabel',AA);
    caxis([0,1]);
    colorbar;
    title('Crest-Trough');
    for i = 1:length(AA),
for j = 1:length(BB),
text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
end;
end;
%set(gca,'FontSize',16)
    subplot(2,2,3);
    imagesc(BB,AA,LV_matrix');
    set(gca,'YTick',AA);
    set(gca,'YTickLabel',AA);
    caxis([0,1]);
    colorbar;
    title('L-value');
 for i = 1:length(AA),
for j = 1:length(BB),
text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
end;
end;
  % set(gca,'FontSize',16)
    subplot(2,2,4);
    imagesc(BB,AA,ED_matrix');
    set(gca,'YTick',AA);
    set(gca,'YTickLabel',AA);
    caxis([0,1]);
    colorbar;
    title('Crest');
    for i = 1:length(AA),
    for j = 1:length(BB),
    text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
    end;
    end;
%set(gca,'FontSize',16)
end;